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Abstract 



Data collected by the Pierre Auger Observatory provide evidence for anisotropy 
in the arrival directions of the cosmic rays with the highest energies, which are 
correlated with the positions of relatively nearby active galactic nuclei (AGN) [1]. 
The correlation has maximum significance for cosmic rays with energy greater than 
^ 6 X 10^^ eV and AGN at a distance less than ^ 75 Mpc. We have confirmed the 
anisotropy at a confidence level of more than 99% through a test with parameters 
specified a priori, using an independent data set. The observed correlation is com- 
patible with the hypothesis that cosmic rays with the highest energies originate from 
extra-galactic sources close enough so that their flux is not significantly attenuated 
by interaction with the cosmic background radiation (the Greisen-Zatsepin-Kuz'min 
effect). The angular scale of the correlation observed is a few degrees, which suggests 
a predominantly light composition unless the magnetic fields are very weak outside 
the thin disk of our galaxy. Our present data do not identify AGN as the sources 
of cosmic rays unambiguously, and other candidate sources which are distributed 
as nearby AGN are not ruled out. We discuss the prospect of unequivocal identifi- 
cation of individual sources of the highest-energy cosmic rays within a few years of 
continued operation of the Pierre Auger Observatory. 



1 Introduction 



The identification of the sources of the cosmic rays with the highest energies 
so far detected has been a great challenge ever since the first event with en- 
ergy around 10^^ eV was reported [2j. If the highest-energy cosmic rays are 
predominantly protons and nuclei, only sources which are less than about 
200 Mpc from Earth could contribute significantly to the observed flux above 
6 X 10^^ eV. Protons with higher energies interact with cosmic microwave 
background photons to produce pions [3,4j, which leads to a significant at- 
tenuation of their flux from more distant sources. The energy of light nuclei 
is damped over an even shorter length scale due to photo-disintegration pro- 
cesses [5,6j. If the relatively nearby sources are not uniformly distributed then 
we expect that the arrival directions of the most energetic cosmic rays should 
be anisotropic, as long as deflections imprinted by intervening magnetic fields 
upon their trajectories are small enough that they point back to their place 
of origin. 

The Pierre Auger Observatory [7], has been operating in Argentina and taking 
data in a stable mode since January 2004. The large exposure of the surface 
detectors (SD), combined with accurate energy and arrival direction measure- 
ments, calibrated and verified from the hybrid operation with fluorescence 
detectors (FD), provides an opportunity to find the clues that could lead to 
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an understanding of the origin of the highest-energy cosmic rays. 

AGN have long been considered possible sites for energetic particle produc- 
tion, where protons and heavier nuclei could be accelerated up to the highest 
energies measured so far pi9] . Windows of a few degrees around each known 
AGN lying within 100 Mpc cover a significant fraction - but not most - of 
the sky. We were therefore motivated to search for an excess, as compared to 
expectations for an isotropic fiux, of cosmic rays with arrival directions close 
to AGN. The angular size of the search window should not be limited to the 
instrumental angular resolution, since correlation could exist on larger scales 
due to magnetic defiections, the precise amount of which is unknown. Arrival 
directions of cosmic rays are reconstructed by the SD array with an angular 
accuracy better than 1° above 10^^ eV [lOj. 

We have recently reported [1] the observation of a correlation between the ar- 
rival directions of the cosmic rays with highest energies measured by the Pierre 
Auger Observatory and the positions of nearby AGN from the 12^^ edition of 
the catalogue of quasars and active nuclei by Veron-Cetty and Veron [11]. In 
this article we provide more details about the methods used to demonstrate 
anisotropy based on this correlation, and further analyse its properties and 
implications. 



2 Evidence for anisotropy and correlation with AGN 

2.1 Data set 

The southern site of the Pierre Auger Observatory [7J is located in Malargiie, 
Argentina, at latitude 35.2° S, longitude 69.5° W, and mean altitude 1400 
meters above sea level. The data set analysed here consists of events recorded 
by the Pierre Auger Observatory from 1 January 2004 to 31 August 2007. 
During this time, the size of the Observatory increased from 154 to 1388 sur- 
face detector stations. We consider events with reconstructed energies above 
40 EeV (1 EeV = 10^^ eV) and zenith angles smaller than 60°. The quality cut 
implemented in the present analysis requires that at least five active nearest 
neighbours surround the station with the highest signal when the event was 
recorded, and that the reconstructed shower core be inside an active equilat- 
eral triangle of detectors. 

The event direction is determined by a fit of the arrival times of the shower 
front at the SD. The precision achieved in the arrival direction depends on the 
clock resolution of each detector and on the fiuctuations in the time of arrival of 
the first particle fT2^ . The angular resolution is defined as the angular aperture 
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around the arrival directions of cosmic rays within which 68% of the showers 
are reconstructed. This resolution has been verified experimentally jl3lll0j . 
Almost all events with energies above 10 EeV trigger at least 6 surface stations 
and have an angular resolution better than 1° p!3llT0] . 

The energy of each event is determined in a two-step procedure. The shower 
size aS, at a reference distance and zenith angle, is calculated from the signal 
detected in each surface station and then converted to energy using a linear 
calibration curve based on the fiuorescence telescope measurements [Hj . The 
uncertainty in S resulting from the adjustment of the shower size, the con- 
version to a reference angle, the fiuctuations from shower-to-shower and the 
calibration curve amounts to about 18%. The absolute energy scale is given by 
the fiuorescence measurements and has a systematic uncertainty of 22% [15]. 
There is an additional uncertainty in the energy scale for the set of high energy 
events used in the present analysis due to the relatively low statistics available 
for calibration in this energy range. 

2. 2 Exposure 

The integrated exposure for the event selection described in the previous sec- 
tion amounts to 9,000 km^ sr yr. Note that analyses involving a fiux calcula- 
tion, such as the estimate of the cosmic-ray spectrum [Hj, use stricter selection 
criteria which would amount to an exposure of about 7000 km^ sr yr for the 
same data period. 

The surface detector array has full acceptance for events with energy above 
3 EeV [16j. Above this energy the detection efiiciency is larger than 99% and it 
is nearly independent of the direction of the shower axis defined by the zenith 
angle (9) with respect to the local vertical and azimuth (0) with respect to 
the South. Thus, above that energy the instantaneous instrument aperture as 
a function of zenith angle is given by : 

A{t) = n{t)aoco8edQdt (1) 

where ag cos 6 is the surface of a unitary cell under the incidence zenith angle 
9 and n{t) is the number of active such cells as a function of time. The number 
n{t) is recorded every second by the trigger system of the Observatory and 
refiects the array growth as well as the dead period of each detector. Such 
recording allows for a precise knowledge of our aperture at any moment in 
time. 

The instrument exposure above a certain energy E may be further affected by 
the conversion of the measured signal at ground to energy (this dependence 
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is not included in Eq. 1 above). For a given energy E the ground signal vary 
depending on the atmospheric conditions (e.g. through the variations of the 
Moliere radius) [17J. If the signal to energy conversion does not correct for 
these small variations, of order a few %, as it is the case in our analysis, the 
aperture above a certain uncorrected energy will depend on the atmospheric 
conditions. 

Over the period from 1 January 2004 to 31 August 2007 the integration of 
the time dependences from the array growth and dead time together with the 
atmospheric variations introduce a modulation of the exposure as a function of 
celestial right ascension (RA) of less than 1%. For the purpose of our analysis, 
where the total number of events considered is less than 100, such modulation 
is negligible and the resulting RA dependence can be safely ignored. 

Hence our exposure only depends on the celestial declination 5 and can be 
derived from the relation sin 5 = cos sin A — sin cos A cos 0, where A is the 
latitude of the Observatory. 



2.3 Search method 

We denote by p the probability that an individual event from an isotropic 
flux has, by chance, an arrival direction closer than some particular angular 
distance ip from any member of a collection of candidate point sources, p is the 
exposure-weighted fraction of the sky accessible to observation by the Pierre 
Auger Observatory which is covered by windows of radius iJj centred on the 
selected sources. 

The probability P that k or more out of a total of events from an isotropic 
flux are correlated by chance with the selected objects at the chosen angular 
scale is given by the cumulative binomial distribution: 



For this analysis we consider the correlation between cosmic rays and AGN in 
the 12*^ edition of the catalogue of quasars and active nuclei by Veron-Cetty 
and Veron [llj (V-C). This catalogue can not be claimed to contain all existing 
AGN, nor to be an unbiased statistical sample of them. It, however, contains 
the results of a thorough survey of all such objects in the literature. This cata- 
logue contains 85,221 quasars, 1,122 BL Lac objects and 21,737 active galaxies. 
Among these objects, 694 have redshift z < 0.024, a value corresponding to 




(2) 
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a distance smaller than approximately 100 MpcQ At distances greater than 
100 Mpc the catalogue becomes increasingly incomplete and inhomogeneous. 
The V-C catalogue is also particularly incomplete around the galactic plane. 
This is not an obstacle to demonstrating the existence of anisotropy but may 
affect our ability to identify the cosmic-ray sources unambiguously. 

We compute the degree of correlation as a function of three parameters: the 
maximum AGN redshift Zjnax^ the maximum angular separation and the 
lower threshold energy for cosmic rays Eth- Our scan in angular separation ip 
is constrained by the angular resolution of the SD [lOj at the low end (we use 
i^min = 1°) and by the increase in the individual probability p at the high end; 
large i/j push the value of p toward unity, rendering searches for correlation 
above isotropic expectations meaningless. Illustrative values are p = 0.27 and 
p = 0.6 for maximum angular distance '0 = 3° and '0 = 6°, respectively, in 
the case of maximum AGN redshift z^ax = 0.024. Our scan in energies is 
motivated by the assumption that the highest-energy cosmic rays are those 
that are least deflected by intervening magnetic flelds, and that they have a 
smaller probability to arrive from very distant sources due to the GZK eflFect 
[Hill] . The scan in energy threshold is carried out starting with the event with 
the highest energy and adding, one by one, events with successively lower 
energy. 

We scan with the method described above to flnd the minimum value of P, 
given in Eq. [2j Note however that Pmin is not the chance probability that the 
observed arrival directions are isotropically distributed. An estimate of the 
chance probability must incorporate the eflPect of the scan performed upon the 
data. To do so, we build simulated sets, each having the same number of events 
as in the data set, drawn from an isotropic flux in proportion to the relative 
exposure of the Observatory. The chance probability is estimated from the 
fraction of simulated isotropic sets that have, anywhere in the parameter space 
and under the same scan, equal or smaller values of Pmin than the minimum 
found in the data [18\. The result can only be considered an estimate of the 
chance probability, since it depends somewhat on the choice of the range for 
the scan parameters, and does not account for the possibility of dilution due 
to diflFerent scan methods or for scans against different sets of astronomical 
objects. 



2.4 Exploratory scan and anisotropy confirmation 



An exploratory search for correlation between cosmic rays and AGN was con- 



ducted according to the method described in Section 2.3 using data collected 



^ For a redshift z small compared to 1, the distance to an object is approximately 
42 Mpc X (z/0.01) for a Hubble constant Hq = 71 km s"^Mpc"^. 
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from 1 January 2004 to 27 May 2006. This search yielded a minimum probabil- 
ity Pmin for the parameter set: Zjnax = 0.018 {Dj^ax = 75 Mpc), Et^ = 56 EeV 
and %l) = 3.1°, with 12 events among 15 correlated with at least one of the 
selected AGN. For this parameter set, the chance correlation is p = 0.21. Only 
3.2 events were expected to correlate by chance if the flux was isotropic. 

Much of the discussion regarding past evidence for possible anisotropy in the 
distribution of the arrival directions of ultra-high energy cosmic rays has been 
centred on the issue of the impact of trial factors on the statistical signiflcance 
of any potential signal. An accurate measure of the statistical signiflcance of 
some previous reports of anisotropy p!9ll20ll21|[22j . could not be achieved due 
to the posterior nature of the analyses involved. It is only from subsequent ob- 
servations that those claims could be quantitatively evaluated |i23ll24ll25ll26ll27j . 

Therefore, to avoid the negative impact of trial factors in a posteriori anisotropy 
searches, the Pierre Auger collaboration decided that any potentially interest- 
ing anisotropy signal should be tested on an independent data set with param- 
eters specifled a priori. This method was described in [28j where a particular 
set of parameters and sources were proposed and subsequently tested on the 
flrst Auger data set [29] . 

The correlation observed in the exploratory scan motivated the construction 
of a speciflc test to reject or accept the isotropy hypothesis with parameters 
specifled a priori on an independent data set, using exactly the same recon- 
struction algorithms, energy calibration and quality cuts for event selection 
as in the exploratory scan. All details of the prescribed test were documented 
and archived in an internal note. 

The test null hypothesis is isotropy and its statistical characteristics are fully 
deflned by the choice of two probabilities known as the type I and type II 
errors. The type I error [a] is the probability of rejecting the null hypothesis 
incorrectly. In our case this is the probability of declaring our independent 
data set anisotropic when it is not. We have chosen a = 1%. The type II error 
(/?) is the probability of accepting the null hypothesis incorrectly. In our case 
this is the probability of declaring the independent data set isotropic when it 
is not. We have chosen /? = 5%. 

The selection and correlation criteria for the events were chosen according 
to the parameter set that minimised the probability in the exploratory scan 
= 3.1°, Zjnax = 0.018, Eth = 56 EeV). Since we could not predict how many 
events would be required to conflrm the results at a statistically signiflcant 
level from the exploratory scan, we adopted a running prescription (with a 
pre-deflned stopping rule) for conducting a sequential analysis with individual 
tests to be applied after the detection of each subsequent event passing our 
selection criteria. 
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Table 1 

Criteria for our running prescription where N corresponds to the total number 
of events observed at any point during the sequential analysis of up to 34 events 
arriving with energy £^ > 56 EeV. kmin is the minimum number of events within the 
angular window ('0 = 3.1°), and a maximum AGN redshift (zj^ax — 0.018) required 
to reject isotropy with at least a 99% confidence level. This prescription applied to 
data collected after 27 May 2006 was satisfied with = 8 and /c = 6 on 25 May 
2007. 
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Fig. 1. Likelihood ratio i? as a function of the number of events observed in the 
prescribed test. The null hypothesis (isotropy) was rejected at the 99% likelihood 
threshold with 10 events. Shaded regions indicate expectations from isotropy at the 
68% and 95% confidence limit. 

If, in the sequence, one of the individual tests is satisfied, we reject the hypoth- 
esis of isotropy with a confidence level of at least (l-a)=99%. The total length 
of the test sequence (34 events) was determined by the requirement of detect- 
ing a minimum correlation power of 60%, as estimated from the statistics of 
the exploratory scan, within our specified f3 of 5%. In Table 1 we list, for a 
given number of events passing our selection criteria A^, the minimum number 
of events in correlation kjnin necessary to reject the null hypothesis (isotropy) 
with a confidence level larger than 99%, accounting for the sequential nature 
and finite length of the test. Note that for some values of (e.g., 5, 7, 11, 
etc.) there exists no value of k^in that can satisfy the threshold probability 
without also having already satisfied the threshold at a lower value of A^. 

The prescribed test was applied to data collected after 27 May 2006, with ex- 
actly the same reconstruction algorithms, energy calibration and quality cuts 
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for event selection as in the exploratory scan. On 25 May 2007, 6 out of 8 
events correlated, thus satisfying the prescription. In the independent data 
set collected up to 31 August 2007 there are 13 events with energy above 
56 EeV, of which 8 have arrival directions closer than 3.1° from the positions 
of AGN less than 75 Mpc away, with 2.7 expected on average if the arrival di- 
rections were isotropic. The probability for this single configuration to happen 
by chance if the fiux were isotropic (Eq. [2]) is P = 1.7 x 10~^. 

Following our search protocol and based on the independent data set alone, we 
reject the hypothesis of isotropy of the arrival directions of the highest-energy 
cosmic rays with at least 99% confidence level. 

An alternative standard technique in sequential analysis could also have been 
used to monitor the evolution of the correlation signal: the sequential likeli- 
hood ratio test [30II3T] . For the sequential test of AGN correlation, the like- 
lihood ratio R is given by the relative binomial probabilities of the isotropic 
(binomial parameter p = 0.21 in our case) and anisotropic (binomial param- 
eter pi > p) cases. Since pi is not known, we integrate over p < pi < 1 to 
obtain the test ratio, R as defined by Wald [3 1,32 J 

a p'lji - Pi^-' dp, 

pk{i-p)^-k+^ ' ^ ' 



The test rejects the isotropy hypothesis whenever R > {1 — f3)/a (95 in our 
case) with the type I error a = 1% and type II error /? = 5% as previously 
specified. 

This ratio test concluded at the rejection of the isotropy hypothesis when 7 
out of the first 10 events of our independent data set correlated with AGN 
locations (see Figure [!} which also shows the subsequent evolution of the 
signal) . 



3 The AGN correlation signal 



Having determined that an anisotropy exists according to an a priori search 
over an independent subset of the Auger data, we now consider results us- 
ing the full data set (1 January 2004 - 31 August 2007) which allows us to 
obtain a more accurate measurement of the correlation signal. This data set, 
constructed using an updated version of our reconstruction algorithm (see ap- 
pendix A), contains 81 events with energy above 40 EeV and zenith angle 



smaller than 60°, which satisfy the quality criteria given in section 2.1 
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3. 1 Maximum correlation parameters 



Using the method described in Section [ZS] apphed to the full data-set, we per- 
formed a scan within the range of parameters 1° < < 8°, < ^max ^ 0.024 
and Eth > 40 EeV. Catalogue-incompleteness prevents reliable exploration of 
higher redshifts. The scan in maximum angular distance is performed in steps 
of 0.1°, and the scan in maximum redshift z^ax is done in steps of 0.001. 



The minimum probability for the hypothesis of isotropic arrival directions 
is found for the parameter set z^ax = 0.017 {Djj^ax ^ 71 Mpc), = 3.2°, 
and Eth = 57 EeV. These results are statistically consistent with the results 
obtained from the earlier exploratory scan. 



With these selected parameters, we find that 20 out of 27 cosmic-ray events 
correlate with at least one of the 442 selected AGN (292 in the field of view 
of the Observatory), while only 5.6 are expected on average to do so if the 
flux were isotropic = 0.21). The respective cumulative binomial probability 
(Eq. [2]) of achieving this level of correlation from an isotropic distribution is 
Pmin = 4.6 X 10~^. The chance probability that the observed correlation arose 
from an isotropic flux is much larger than Pmim as already discussed in section 



2.3, because a scan was performed over a large parameter space to flnd the 



minimum of P. 



To account for the effects of the scan we built simulated sets each with equal 
number of arrival directions (81 in our case) drawn from an isotropic flux 
in proportion to the relative exposure of the Observatory, and counted the 
fraction of simulated sets which had, anywhere in the parameter space and 
under the same scan, equal or smaller values of Pmin than the minimum found 
in the data [18j. With this procedure, we obtained smaller or equal values of 
Pmin in 10~^ of the simulated sets. 



In Figure [2] we present a sky map, in galactic coordinates, with circles of 
radius 3.2° around each of the arrival directions of the 27 events with energy 
E > 57 EeV detected by the Pierre Auger Observatory, along with asterisks 
at the positions of the 442 AGN with redshift z < 0.017 in the V-C catalogue. 
Each coloured band represents an equal integrated exposure which varies by 
about a factor of 3 between the lightest and darkest band. The number of 
AGN in each of those 6 bands is given in Table [2] The energies and arrival 
directions of the events are listed in Appendix [Aj 
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Table 2 

Number of AGN with z <= 0.017 in each of the exposure bands indicated in Fig. [2j 
Each of the top 6 bands represent 1/6 of the total exposure, the corresponding 
fraction of the whole sky is also indicated. The last declination band represents the 
part of the sky outside the field of view of Auger for zenith angles < 60°. 




Fig. 2. Aitoff projection of the celestial sphere in galactic coordinates with circles 
of 3.2° centred at the arrival directions of 27 cosmic rays detected by the Pierre 
Auger Observatory with reconstructed energies E > 57 EeV. The positions of the 
442 AGN (292 within the field of view of the Observatory) with redshift z < 0.017 
{D < 71 Mpc) from the 12^^ edition of the catalogue of quasars and active nuclei 
[TT ] are indicated by asterisks. The solid line draws the border of the field of view 
for the southern site of the Observatory (with zenith angles smaller than 60°). The 
dashed line is, for reference, the super-galactic plane. Darker colour indicates larger 
relative exposure. Each coloured band has equal integrated exposure. Centaurus A, 
one of the closest AGN, is marked in white. 
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Fig. 3. Probability for the null hypothesis (isotropic distribution) vs. maximum 
angular distance '0 (left), maximum AGN redshift Zjnax (centre), and threshold 
cosmic-ray energy Ef^ (right). In each case the other two parameters are held fixed 
at the values that lead to the absolute minimum probability ('0 — 3.2°, Zfjiax — 

0.017, 

Eth = 57 EeV). 



3.2 Properties of the correlation signal 



In Figure [3] we show one-dimensional plots of the probability P as a function 
of each of the scan parameters with the other two held fixed at the values 
which lead to the absolute minimum probability. 

We note that the energy threshold at which the correlation with nearby AGN 
is maximised, i.e., Efh = 57 EeV, matches the energy range at which the 
flux measured by the Pierre Auger Observatory is ^ 50% lower than would 
be expected from a power law extrapolation of the spectrum measured at 
lower energies [14j. This feature adds support to the interpretation that the 
correlation with relatively nearby sources is evidence for the GZK effect [3,4J, 
as will be discussed in Section 14.41 

Relatively small values of P occur for the energy threshold Eth ^ 60 EeV for 
a range of maximum distances to AGN between 50 Mpc and 100 Mpc and 
for angular separations up to 6°. For instance, there is a local minimum with 
a value P = 8 x 10~^ very close to that of the absolute minimum {Pmin = 
4.6 X 10"^) for the set of parameters i/j = 4.8°, Zmax = 0.013 {D^ax = 55 Mpc). 
With this set of parameters there are 22 events among the 27 with E > 57 EeV 
that correlate with at least one of the 310 selected AGN, while only 7.4 were 
expected, on average, to do so by chance if the flux was isotropic {p = 0.28). 
With limited statistics, the parameters that minimise the probability P should 
only be taken as indicative values of the relevant correlation scales. 

AGN catalogues are likely to be incomplete near the galactic plane, where 
extinction from dust in the Milky Way reduces the sensitivity of observations. 
Moreover, cosmic rays that arrive close to the galactic plane are likely to 
have been deflected by the magnetic field in the disk more than those which 
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Fig. 4. Distribution of angular separations to the closest AGN within 71 Mpc. The 
6 events with |6| < 12° have been shaded in grey. The average expectation for an 
isotropic flux is shown as the dashed line histogram. 

arrive with higher galactic latitudes. These effects could have some impact 
upon the estimate of the strength and of the parameters that characterise the 
correlation. Catalogue incompleteness would weaken the measured strength of 
a true correlation. 

In Figure |4] we plot the distribution of angular separations between the arrival 
directions of the 27 highest-energy events and the position of the closest AGN 
with redshift z < 0.017. On this graph the 6 events with galactic latitudes 
\b\ < 12° have been shaded in grey. The two distributions are clearly distinct, 
a likely consequence of the incompleteness of the V-C catalogue at low galactic 
latitudes. The dashed line is, for comparison, the distribution expected, on 
average, from an isotropic flux modulated by the relative exposure of the 
Observatory. 

We have performed a scan limited to events with galactic latitudes \b\ > 12°. 
The minimum probability for the hypothesis of isotropic arrival directions 
occurs for the same parameters as without the cut in the galactic plane {ip = 
3.2°,z^ax = O.On^Eth = 57 EeV). The cut increases the strength of the 
correlation {Pmin = 1-7 x 10~^°). 19 out of 21 arrival directions correlate with 
AGN positions while 5.0 are expected to do so by chance if the flux were 
isotropic. In other words, 5 of the 7 events which do not correlate with AGN 
positions arrive with galactic latitudes \b\ < 12°. 

A distribution of arrival directions of cosmic rays that shows evidence of 
anisotropy by correlation with a set of astrophysical objects is also expected 
to show evidence for anisotropy by auto-correlation. The degree of auto- 
correlation in the set of the 27 events with E > 57 EeV is shown in Figure [5| 
where we plotted the number of event pairs with angular separation smaller 
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Fig. 5. Number of pairs as a function of maximum separation angle a for the 27 
events with E > 57 EeV (points) and average expectation for an isotropic flux. 
The error bars on the isotropic expectations represent the 90% confidence limit 
dispersion. 



than a given value. Points represent the number of pairs in the data. Also 
shown are the mean number of pairs expected in simulated isotropic sets of 
27 directions, distributed in proportion to the exposure of the Observatory. 
The error bars represent the dispersion of 90% of the simulations. Significant 
departures from isotropy are seen to occur at intermediate angular scales, be- 
tween 9° and 22°. This may be the consequence of a combination of clustering 
of events from individual sources in addition to effects of the non-uniform 
distribution of the sources themselves [26j . 

To compare the auto-correlation function of the data to that expected from 
the AGN distribution in the V-C catalogue we must restrict ourselves to the 
regions where the catalogue is reasonably complete, e.g., outside of the galactic 
plane. In Figure |6] we plot the number of pairs in the data as a function of 
the separation angle restricted to the 21 events with E > 57 EeV and galactic 
latitudes \b\ > 12°. Also shown is the average distribution expected in sets of 
21 directions chosen at random (in proportion to the relative exposure of the 
Observatory) from the positions of AGN in the V-C catalogue with redshift 
z < 0.017 and \b\ > 12°. The error bars in the plots indicate the results in 
90% of the simulated sets. The distribution of pairs in the data are in all cases 
within those results. 
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Fig. 6. Number of pairs as a function of maximum separation angle a for the subset 
of 21 events with E > 57 EeV and |6| > 12° (points) and average expectation for 
AGN in the V-C catalogue with z < 0.017 and the same cut in galactic latitude. The 
error bars on the AGN expectations represent the 90% confidence limit dispersion. 

Anisotropy in the distribution of arrival directions can also be characterised 
by the multipole moments of an expansion in spherical harmonics l^^rn- With 
partial coverage of the sky the application of this method is not straightfor- 
ward [33] . A method has been developed [34] to cope with partial sky coverage 
that allows a bound to be placed, with a given confidence level, to the mini- 
mum multipole order necessary to account for the observed distribution. The 
application of this method to the arrival directions of the events with energies 
above 50 EeV in our data set shows that a multipole order ^ > 5 is necessary 
to describe their distribution with 99% confidence level. The arrival directions 
of events with energy between 40 EeV and 50 EeV do not show evidence of 
anisotropy with this method. This result is consistent with the sharp decrease 
of the correlation with AGN positions at lower energies. 



4 Discussion 



In the previous sections we have demonstrated the anisotropy of the highest- 
energy cosmic rays and have derived a set of parameters that maximises the 
correlation with the AGN sample from the 12th edition of the V-G catalogue. 
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In this section, based on those parameters and the data set that maximises 
the correlation, we discuss some possible implications of the observed signal. 

4.1 Lower limit on the number of sources 

If ultra high energy cosmic rays come from a large number of dim sources, 
the number of pairs of events (doublets) coming from one source is expected 
to be much smaller than the number of singlets. On the other hand, if they 
come from a small number of bright sources, the ratio of doublets to singlets 
is expected to be larger. It is then possible to put a lower limit on the number 
of sources based on the ratio of doublets to singlets. The minimum number 
of sources, S', results for the case in which all the sources have the same 
apparent luminosity [35j. If sources are steady, cosmic rays accelerated by 
one source at different times are statistically independent and the detection 
can be considered as a Poisson process. Then the probability that one source 
accelerates n particles is given by P{n) = v^e~^ /n\^ with v the mean number 
of events expected from one source. The mean number of expected singlets 
from S sources is ni = x P(l) = Siye~^ and that of doublets is n2 = 
S X P{2) = Sv^e~^ /2. It is thus possible to estimate the number of sources S 
as: 

S ^ ^e(2-2)M (4) 

2712 ^ ^ 

Within the 27 highest-energy events there are 6 pairs with separation smaller 
than the correlation angular scale of 6°, while 1.6 are expected by chance in 
an isotropic flux. Taking 77,2 = 6 — 1.6 = 4.4 and ni = 27 — 2n2 = 18.2, we 
obtain a lower limit for the number of sources S' > 61. 

Note that this is a bound for mean expectations, but could have large fluc- 
tuations with the present small statistics. Also, it was derived under the un- 
realistic assumption of equal flux on Earth for all sources. Assuming instead 
equal intrinsic luminosity in cosmic rays the mean number of sources becomes 
larger by a factor of order ni/n2 [35j. The lower bound could also increase if 
the sources had signiflcant clustering of their own on the same angular scale 
as the clustering of events. In either case, this lower limit does not contradict 
the hypothesis that nearby AGN are the sources. 

4-2 Signal dependence on energy 

We have studied the dependence on energy of the correlation of our data 
set with the AGN from the V-C catalogue. In an approach similar to the 
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one developed in [36] we constructed a smoothed density map from the V-C 
catalogue and used it to compute the log-likelihood of any event sample. We 
then compared the result for the data to that for simulated samples of the 
same size, either drawn from an isotropic distribution of arrival directions or 
from the smoothed density map itself, in both cases modulated by the relative 
acceptance of the Observatory. 

From the density map of the V-C catalogue, smoothed on a given angular 
scale and limited to within 100 Mpc, it is possible to calculate the average 
log-likelihood of an event sample : 

1 ^ 

LL=M^ ^""SiPk) (5) 

k=i 

where the sum runs over the events of the sample under consideration and 
pk is the map density in the direction of the event k. 

In Figure [T] we present the values of LL (dots) from our data in 10 EeV energy 
intervals. We used events with galactic latitudes |6| > 12°, and compared 
them with the average expected from samples of similar size and generated 
either according to an isotropic distribution (bottom (blue) line) or to the 
distribution of AGN within 100 Mpc smoothed with Gaussian windows of 
2° (top (red) line). The dashed lines represent the 1 a interval around the 
mean value, based on the statistics of the real data. The data are compatible 
with a distribution that follows that of AGN at high energies with an abrupt 
transition towards an isotropic distribution below 60 EeV. 

4^3 Source identification 

The high degree of correlation that we observe can certainly serve as a strong 
indication that AGN may well be the sources of ultra-high energy cosmic rays. 
However this result is not yet a proof. 

In particular, we know that the distribution of matter (visible and dark) in 
our local Universe (within 100 Mpc) is strongly non-uniform, and that AGN 
are correlated with this non-uniformity. We are therefore motivated to see if 
our data provide sufficient information to determine if the correlation signal 
is unambiguously associated with AGN or if they are simply acting as tracers 
of some other set of sources with a similar distribution. 

In a specific study we have compared the arrival directions of our highest 
energy events with the distribution expected from different source models 
using a likelihood test. The test compared isotropic distributions, distributions 
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Fig. 7. Comparison of the average log-likelihood per event in different energy in- 
tervals (calculated using the smoothed distribution of AGN) between real data and 
simulated samples (bottom [blue line] isotropic distribution, top [red line] AGN 
distribution from V-C catalog). 

of galaxies from the IRAS PSCz [37j catalogue and distributions of AGN from 
the the V-C catalogue, at different angular scales and using different horizon 
depths (see 4.4) . In addition to confirming that our data are incompatible with 
an isotropic distribution these studies showed that they are best represented 
by a small angular smoothing within a relatively nearby horizon (100 Mpc) of 
the AGN from the V-C catalogue [36] . 

The local spatial distribution of AGN is correlated with the distribution of 
other astronomical objects that are potential sources of cosmic rays, such as 
rich clusters of galaxies and star-burst galaxies (which could host a large frac- 
tion of gamma ray bursts). We have tested, with the same scan method that 
we applied to the V-C catalogue, whether there is a significant correlation with 
the Abell clusters of galaxies (an optical survey) [38] and with the clusters in 
the X ray surveys REFLEX [39j and NORAS |ffi]. We also tested for correla- 
tion against a catalogue of starburst galaxies [^T^ . We did not find significant 
correlation in these searches with our present data set. 



It is plausible that only a subclass of AGN in the V-C catalogue are the sources 
of the highest-energy cosmic rays. With our present relatively small data set 
it is difficult to pinpoint distinctive properties of the AGN that are close to 
their arrival directions, or to draw firm conclusions about patterns in their 
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redshift distribution. It is worth noting, as is clearly visible in Figure [2} the 
striking alignment of several events close to the super-galactic plane. Two of 
the events have arrival directions less than 3° away from Centaurus A, one of 
the closest AGN. 



44 The GZK horizon 

The correlation observed is consistent with the hypothesis that the highest- 
energy cosmic rays that arrive on Earth are predominantly produced in rel- 
atively nearby AGN, within the distance over which the GZK effect j3ll4ll42] 
does not significantly attenuate their fiux. 

The "GZK horizon" may be defined as the distance from the Earth which 
contains the sources that produce 90% of the protons that arrive with ener- 
gies above a given threshold. Under the idealisation of uniformly distributed 
sources of equal intrinsic cosmic-ray luminosity and a conventional spectral 
index, the horizon computed in the continuous energy loss approximation is 
about 90 Mpc for protons that arrive with energies above 80 EeV and about 
200 Mpc for energies above 60 EeV [43j . Deviations of the horizon scale from 
the estimates above are expected, in particular due to local departures of 
the sources from uniformity in spatial distribution, intrinsic luminosity, and 
spectral features. 

The largest departure from isotropic expectations (minimum value of the prob- 
ability P) in the complete data set was found to be due to correlation with 
AGN at a distance smaller than 71 Mpc and for cosmic rays with energies 
above 57 EeV. However, relatively small values of P occur for this energy 
threshold for a range of maximum distances to AGN between 50 Mpc and 
100 Mpc. 

If these numbers were to be taken at face value, an upward shift in the energy 
calibration of ^ 30%, as suggested in some simulations of the reconstruc- 
tion of the shower energies [44j, would lead to a better agreement between 
the maximum AGN distance Dmax that minimises the probability P and the 
theoretical expectations based on the idealised GZK attenuation. However, 
while we expect i^rnax to be comparable to the GZK horizon scale, the re- 
lation is not a simple one. An accidental correlation with foreground AGN 
different from the actual source may induce some bias in the value of i^rnax 
toward smaller maximum source distances. The P minimisation method has 
non-uniform sensitivity over the range of parameters explored. Incompleteness 
of the V-C catalogue prevents a reliable exploration beyond approximately 
100 Mpc. Furthermore, as mentioned above, a realistic estimate of the hori- 
zon scale depends on several unknown features. A large local over-density 
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of sources would reduce its value. The distribution of intensity and spectral 
features of the dominant sources also has an effect on the horizon scale. 

Regarding the possibility that the cosmic rays injected at the sources are 
heavy nuclei, attenuated mainly by photo-disintegration processes, one may 
note that nuclei of the iron group have horizons only slightly smaller than the 
proton horizons, but intermediate mass nuclei {A 20-40) have significantly 
smaller horizons (e.g., the horizon for a threshold energy of 60 EeV is about 
60 Mpc for ^^Si nuclei [43]). The smaller horizon for decreasing nuclear mass 
is due to the corresponding decrease in the threshold required to excite the 
giant-dipole resonance for photo-disintegration. 



4^5 Effect of the magnetic fields 

A cosmic ray with charge Ze that travels a distance in a regular magnetic 
field B is defiected by an angle S given by 



If the regular galactic magnetic field has a strength of a few /xG with a coher- 
ence scale of order ^ 1 kpc, as in some models [45j, the defiection is expected 
to be a few degrees for protons with > 60 EeV. In such models, the angular 
scale of the correlation we observed is consistent with the size of the defiections 
expected to be imprinted upon protons by the galactic magnetic field. 

The precise amount of the defiection is very dependent on each specific arrival 
direction. We have evaluated numerically the defiections imparted in a conven- 
tional regular galactic magnetic field model |45j for sets of arrival directions 
uniformly distributed according to the Pierre Auger Observatory relative ex- 
posure. Anti-particles were backtracked in the galactic regular magnetic field 
to a distance of 20 kpc away from the Galactic Centre (where the field strength 
is already very small). At this point the angle between the initial (as measured 
on the Earth) and final velocity vectors was calculated. The result is shown in 
the left panel of Figure [s] in the special case of the BSS-S modep] (without B, 
component) for = 60 EeV protons. The defiections scale approximately as 
Z X (60 EeV/E) for other energies and electric charges (the scaling is rigorous 
only for very small defiections and a uniform field). In the right panel of the 
figure we show the distribution of defiections for protons in the case of the 27 



^ In fact, we have smoothed the original BSS-S model of [45j as described in [|6] 
in order to avoid the discontinuities present in the original model. 




(6) 
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angle (degrees) 



angle (degrees) 



Fig. 8. Distribution of the deflections for protons in the BSS-S model of the galactic 
magnetic field. Left panel: 1000 directions drawn from an isotropic flux in proportion 
to the exposure of the Observatory, for £" = 60 EeV. Right panel: deflections of the 
27 arrival directions of the observed events with E >57 EeV . 

arrival directions of the events with E > 57 EeV, as computed for each using 
its reconstructed energy. 



Models of the regular component of the galactic magnetic field outline 
its basic features, but cannot be expected to provide a complete picture nor a 
realistic value for every direction. It is, for example, possible to do the exercise 
of "correcting" the observed arrival directions to undo the deviation imparted 
by the galactic magnetic field, but current models are not expected to be ac- 
curate enough to allow us to draw reliable conclusions from such analyses. 
Nonetheless, the results shown in Figure |8] provide a reasonable estimate of 
the typical deflections to be expected. They are consistent with the angular 
scale of the observed correlation with AGN. Therefore, if the BSS-S model 
is a fair representation of the general features of the regular galactic mag- 
netic field, then the correlation observed in the data would be unlikely if the 
primary composition of the cosmic rays reaching us were much heavier than 
protons. Note that this does not preclude the possibility that the source emits 
heavy nuclei, which could disintegrate along their journey, so that the lighter 
fragments are those deflected by the galactic magnetic field. 

It will be possible in the future to compare various models for the galactic mag- 
netic field with the pattern of orientation and size of the deviation between the 
observed arrival directions of the events and potential AGN sources. However, 
our present data set is not large enough to perform such an analysis reliably. 

The angular scale of the observed correlation also implies that intergalactic 
magnetic fields along the line of sight to the sources do not in general deviate 
cosmic-ray trajectories by much more than a few degrees. The root-mean- 
square deflection imprinted upon the trajectories of cosmic rays with charge 
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Ze as they travel a distance in a turbulent magnetic field with coherence 
length Lc is 



^ eOEeV B,^, I D I 
0,^,-4 l0-9GVl00MpcVlMpc 



There is no measurement of the intergalactic magnetic fields except at the 
centres of rich galaxy clusters. Numerical simulations of those fields give a 
wide range of possible deflections from negligible [48j to very large [49]. The 
correlation observed can be used to constrain models of turbulent intergalac- 
tic magnetic fields, which must be such that in most directions Brms\fL'c ^ 
lO-^G^M^ within our "GZK horizon". 

Finally, there are possible biases in the determination of the relevant angular 
scale of the deflection. The active galaxy closest to the arrival direction of a 
cosmic ray is not necessarily the source responsible for it. This could lead to 
an underestimate of the deflection involved. In principle it could also happen, 
due to catalogue incompleteness or because the source is something else, that 
an AGN found to correlate with a particular event is further away than the 
actual source, thus overestimating the deflection angle. 



J1..6 The acceleration sites 



Acceleration sites in the active galaxies that correlate with events above 
57 EeV are promising candidate sources of high energy cosmic rays, but other 
possible sites cannot be ruled out with the present limited statistics. (For a 
recent summary of proposed acceleration sites see, e.g., [50].) The observed 
correlation shows that ultra-high energy cosmic ray sources are extra-galactic 
with an angular distribution similar to that of AGN within ^ 71 Mpc and that 
the primaries are most likely protons that suffer losses due to interactions with 
the cosmic background radiation. These results rule out models for the origin 
of cosmic rays that place observed sources predominantly in our Galaxy, such 
as galactic compact objects (young neutron stars [51j, pulsars [52], and black 
holes), and gamma-ray bursts [53]. Models where sources are located in the 
galactic halo are also ruled out, such as the decay of super-heavy dark matter 
particles [5Hl55ll56j . which are already highly constrained by the Auger limit 
on the fraction of photon primaries at high energies [57,58j. Top-down models 
based on topological defects [59ll60ll61j need to have a spatial distribution con- 
sistent with the local matter distribution to avoid being excluded as significant 
sources of cosmic rays. Such models are also constrained by the photon and 
neutrino limits [58j. 

The large-scale structure distribution of matter, which is traced by normal 
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galaxies, has a similar spatial distribution to the local AGN. Therefore, accel- 
eration sites in galaxies with inactive nuclei cannot be excluded at the present, 
including those based on extra-galactic compact objects [62j, quasar remnants 
[63] . galactic winds in star-bursts galaxies [64j, and gamma-ray bursts [65,66]. 
In contrast, acceleration models in massive clusters of galaxies, such as cluster 
accretion shocks [671168] . are challenged by the observed correlation. Massive 
clusters are rare within 100 Mpc when compared to the number of observed 
events and there is a paucity of events from the direction of Virgo, the nearest 
sizeable cluster of galaxies. 



AGN have long been suggested as likely accelerators of cosmic rays [8ll9ll69] . 
The case for active galaxies as likely sources is based on the power available 
from the central black hole. AGN are powered by the accretion of matter 
onto a super-massive black hole (with masses in the range 10^ to 10^ Mq) 
at the centre of the galaxy. A number of mechanisms have been proposed 
that utilise different regions and properties of this system to accelerate cosmic 
rays to ultra-high energies. Acceleration based on the central regions face the 
challenge of energy losses in the radiation field that surrounds the central black 
hole and accretion disk. Alternative acceleration sites include jets [ 701I7T] and 
radio lobes [72] that are associated with the most luminous AGN. 



AGN with prominent radio lobes are rare and do not follow the observed 
spatial distribution of the observed correlated AGN. The one exception is 
Gen A, at only 3.4 Mpc [73], which has been proposed as a site for cosmic-ray 
acceleration [74j. It displays jets, radio lobes which extend over a scale of about 
10° along the super-galactic plane, and a variable compact radio nucleus. Two 
events correlate with the nucleus position while several lie in the vicinity of the 
radio lobe extension along the super galactic plane (see Figure [2]) . The most 
prominent radio galaxy in our GZK neighbourhood is M87, which does not 
correlate with any observed event above 57 EeV thus far but the coverage of 
the southern Auger Observatory is almost a factor of 3 lower in this direction 
of the sky than in the direction of Gen A. Of the remaining 18 correlating 
events, 15 have Seyfert galaxies as the closest AGN in angular separation. 



A significant increase in ultra-high energy cosmic-ray statistics combined with 
searches for counterparts in a multi- wavelength and multi-messenger campaign 
should improve our ability to distinguish if AGN are the sources of cosmic 
rays or tracers of the sources. If future data select AGN as hosts of cosmic-ray 
accelerators, the type of AGN selected, together with spectral and composi- 
tion information, should help distinguish between proposed AGN acceleration 
mechanisms. 
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5 Conclusions 



Anisotropy has been established with more than 99% confidence level in the 
arrival directions of events with energy above ^ 60 EeV detected by the Pierre 
Auger Observatory. These events correlate over angular scales of less than 6° 
with the directions towards nearby (D < 100 Mpc) AGN. 

The observed correlation demonstrates the extra-galactic origin of the highest- 
energy cosmic rays. It is consistent with the hypothesis that cosmic rays with 
energies above ^ 60 EeV are predominantly protons that come from AGN 
within our "GZK horizon". This provides evidence that the observed steep- 
ening of the cosmic-ray spectrum at the highest energies is due to the "GZK 
effect" , and not to acceleration limits at the sources. 

It is possible that the sources are other than AGN, as long as their local 
distribution is sufficiently correlated with them. Unequivocal identification of 
the sources requires a larger data set, such as the Pierre Auger Observatory 
will gather in a few years. In particular, one could use the fact that angular 
departures of the events from an individual source due to magnetic defiections 
should decrease in inverse proportion to the energy of the cosmic ray. The 
observation of such angle/energy correlation in clusters of events could be 
exploited to locate the source position unambiguously with high accuracy. 
This could also provide at the same time valuable and unique information 
about the magnetic fields along the line of sight. 

We have shown that astronomy of charged particles is indeed feasible at the 
highest energies and that in the next few years we can hope for unambiguous 
identification of sources of cosmic rays. 
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A Event List 



Here we list the 27 events recorded from 1 January 2004 until 31 August 2007 
with energy in excess of 57 EeV. We have indicated the date of observation 
(year and Julian day), the zenith angle, the shower size at 1000 m from the 
core S'(IOOO), the energy in EeV, the equatorial coordinates (RA, Dec) and 
the galactic coordinates (Longitude, Latitude). Events that correlate within 
3.2° of AGN with redshift z < 0.017 are marked with a star. The dashed 
horizontal line indicates the beginning of the prescribed test of section 2. 

The quoted energy is derived from a calibration procedure where the shower 
size is compared to the energy measured by the FD. This energy calibration, 
based on the sample of hybrid events analysed at the time of the prescription, 
was used for the whole data set for consistency. The smaller uncertainty on 
the energy calibration curve expected from the increased statistics of hybrid 
events, as well as improvements in the systematic uncertainty of the FD en- 
ergy scale, may lead to revised energies in future publications. Thus, we also 
include the shower size S(IOOO) at 1000 m from the reconstructed core, as 
it is the shower parameter that is directly measured from the individual SD 
signals in the event. This parameter is almost independent of the shower lat- 
eral distribution function used in the reconstruction procedure (within 10%). 
The uncertainty in S resulting from the adjustment of the shower size, the 
conversion to a reference angle, the fluctuation from shower to shower and the 
calibration curve amounts to about 18%. The absolute energy scale is given by 
the fluorescence measurements and has a systematic uncertainty of 22% [T5]. 

Note that the energies and arrival directions given in this list correspond to 
the analysis of the full data set which used a slightly different reconstruction 
package than the one used for the original scan and the prescribed test of 
section 2. In particular improvements made in the SD tank calibration have 
very slightly modifled the energy and arrival directions. 

If one were to apply the prescribed parameters to this particular reconstruc- 
tion, the prescription would have been fulflUed earlier with the event 2007-069 
(10 March 2007, marked with a double star in the table), with 5 events in 
correlation out of 6 above 56 EeV. 
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Year Julian day 9 5(1000) E (EeV) RA Dec Longitude Latitude 



2004 


125 


47.7 


252 


70 


267.1° 


-11.4° 


15.4° 


8.4° 


2004 


142 


59.2 


212 


84 


199.7° 


-34.9° 


-50.8° 


27.6° * 


2004 


282 


26.5 


328 


66 


208.0° 


-60.3° 


-49.6° 




2004 


339 


44.7 


316 


83 


268.5° 


-61.0° 


-27.7° 


-17.0° * 


2004 


343 


23.4 


323 


63 


224.5° 


-44.2° 


-34.4° 


13.0° * 


2005 


54 


35.0 


373 


84 


17.4° 


-37.9° 


-75.6° 


-78.6° * 


2005 


63 


54.5 


214 


71 


331.2° 


-1.2° 


58.8° 


-42.4° 


2005 


81 


17.2 


308 


58 


199.1° 


-48.6° 


-52.8° 


14.1° * 


2005 


295 


15.4 


311 


57 


332.9° 


-38.2° 


4.2° 


-54.9° * 


2005 


306 


40.1 


248 


59 


315.3° 


-0.3° 


48.8° 


-28.7° * 


2005 


306 


14.2 


445 


84 


114.6° 


-43.1° 


-103.7° 


-10.3° 


2006 


35 


30.8 


398 


85 


53.6° 


-7.8° 


-165.9° 


-46.9° * 


2006 


55 


37.9 


255 


59 


267.7° 


-60.7° 


-27.6° 


-16.5° * 


2006 


81 


34.0 


357 


79 


201.1° 


-55.3° 


-52.3° 


7.3° 


2006 


185 


59.1 


211 


83 


350.0° 


9.6° 


88.8° 


-47.1° * 


2006 


296 


54.0 


208 


69 


52.8° 


-4.5° 


-170.6° 


-45.7° * 


2006 


299 


26.0 


344 


69 


200.9° 


-45.3° 


-51.2° 


17.2° * 


2007 


13 


14.3 


762 


148 


192.7° 


-21.0° 


-57.2° 


41.8° 


2007 


51 


39.2 


247 


58 


331.7° 


2.9° 


63.5° 


-40.2° * 


2007 


69 


30.4 


332 


70 


200.2° 


-43.4° 


-51.4° 


19.2° ** 


2007 


84 


17.3 


340 


64 


143.2° 


-18.3° 


-109.4° 


23.8° * 


2007 


145 


23.9 


392 


78 


47.7° 


-12.8° 


-163.8° 


-54.4° * 


2007 


186 


44.8 


248 


64 


219.3° 


-53.8° 


-41.7° 


5.9° 


2007 


193 


18.0 


469 


90 


325.5° 


-33.5° 


12.1° 


-49.0° * 


2007 


221 


35.3 


318 


71 


212.7° 


-3.3° 


-21.8° 


54.1° * 


2007 


234 


33.2 


365 


80 


185.4° 


-27.9° 


-65.1° 


34.5° 


2007 


235 


42.6 


276 


69 


105.9° 


-22.9° 


-125.2° 


-7.7° 
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